Population differentiation and intraspecific genetic admixture in two Eucryptorrhynchus weevils (Coleoptera: Curculionidae) across northern China

Abstract Increasing damage of pests in agriculture and forestry can arise both as a consequence of changes in local species and through the introduction of alien species. In this study, we used population genetics approaches to examine population processes of two pests of the tree‐of‐heaven trunk weevil (TTW), Eucryptorrhynchus brandti (Harold) and the tree‐of‐heaven root weevil (TRW), E. scrobiculatus (Motschulsky) on the tree‐of‐heaven across their native range of China. We analyzed the population genetics of the two weevils based on ten highly polymorphic microsatellite markers. Population genetic diversity analysis showed strong population differentiation among populations of each species, with F ST ranges from 0.0197 to 0.6650 and from −0.0724 to 0.6845, respectively. Populations from the same geographic areas can be divided into different genetic clusters, and the same genetic cluster contained populations from different geographic populations, pointing to dispersal of the weevils possibly being human‐mediated. Redundancy analysis showed that the independent effects of environment and geography could account for 93.94% and 29.70% of the explained genetic variance in TTW, and 41.90% and 55.73% of the explained genetic variance in TRW, respectively, indicating possible impacts of local climates on population genetic differentiation. Our study helps to uncover population genetic processes of these local pest species with relevance to control methods.


| INTRODUC TI ON
In recent years, many new pests have damaged agricultural production and forestry operations (Riegler, 2018). These pests can be local species or introduced species from other distribution ranges. The local pests show increases in damage due to changes in the environment and management procedures or genetic adaption to local environments (Hoffmann & Sgro, 2011;Riegler, 2018); these include the direct effects of climate changes and management responses to them as well as the evolution of pesticide resistance (Bergé & Ricroch, 2010;Elbert et al., 2019). The spread of species is often facilitated by human activities (Bradshaw et al., 2016;Simberloff et al., 2013). In the past decade, human trade has increasingly promoted the movement of species beyond their historical distribution (Campagnaro et al., 2018;Gippet et al., 2019;Liebhold & Tobin, 2008). Identifying the possible sources of a new pest can provide essential information on developing control methods (Kirk et al., 2013a).
Population genetics approaches can help to reveal the population history of species and then infer whether the species is local or introduced outside of their native range (Antoine et al., 2017;Boissin et al., 2012b;Fraimout et al., 2017). For a pest species in its native range, genetic differentiation of populations is usually aligned to geographic distance (Wei et al., 2015), although the same pattern can be found in introduced species with stepping-stone expansion (Cao et al., 2016;Lundhagen et al., 2013). For a newly introduced species, there is usually a lack of population differentiation or genetic isolation by geographic distance (Cao et al., 2017). The introduced populations usually share the same genetic cluster with their source populations (Beichman et al., 2018;Kirk et al., 2013b).

Multiple introductions in different source populations can all lead
to structured populations in the introduced area (Cao et al., 2017).
Thus, based on patterns of population genetic differentiation and the connectivity among populations, it is possible to infer the population history of a species and provide insights into our understanding of newly outbreaking pests.
The tree-of-heaven trunk weevil (TTW), Eucryptorrhynchus brandti (Harold) and the tree-of-heaven root weevil (TRW), E. scrobiculatus (Motschulsky), are the only two wood-boring species of Eucryptorrhynchus in China (Liu, 2016). The two species both only feed on the tree-of-heaven, Ailanthus altissima (Mill.) Swingle and its varieties A. altissima Qiantouchun in China and cause serious damage. They spend their larvae stage in the trunk (TTW) or root (TRW) of the host tree. There is a strict reproductive isolation between these two species of weevil. TRW lay eggs in the soil near the roots of Ailanthus altissima, while TTW lay eggs in the trunk of Ailanthus altissima (Zhang et al., 2019). In the past decades, TTW and TRW have caused increasing damage to A. altissima and become important pests in northern China (Wu et al., 2016) especially the Ningxia Hui Autonomous Region. This phenomenon is coincident with the widespread use of A. altissima as an ornamental tree planted along the sides of roads in northern China. In North America, TTW was used as a potential biological control agent to control A. altissima where it is an invasive plant (Herrick et al., 2012). Identifying the population evolutionary processes of these weevils may provide insights into local and more widespread movement of the two weevil pests, helping the development of control methods where these weevils are considered pests and ways of promoting their dispersal where the weevils are considered as biological control agents.
In this study, we examined the population genetic diversity and structure of the two weevils to infer their population history.
Comparative study of the two weevils allows us to check if the two co-occurring pests show similar population genetic processes. We hypothesized that the increasing damage by the two species in China is related to human-mediated dispersal of seedings, based on the low flight ability of the two species and their life cycle. We also tested whether local climate can have an impact on the population genetic structure of the two species.

| Sample collection and genotyping
We collected 13 TTW and 11 TRW populations across their native range of China (Table 1, Figure 1). The two weevils are the primary pests of the tree-of-heaven, and we sampled sites separated by variable distances (30-2100 km). To avoid the collection of siblings, one adult individual was collected from each tree. The samples were kept in 100% ethanol and stored at −80°C before DNA extraction. Total genomic DNA was extracted from the legs of each adult using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions.
Ten microsatellite loci from Zhang et al. (2021) were used in our study (Table S1) and genotyped using the same methods. The size of amplified PCR products was determined using an ABI 3730xl DNA Analyzer with GeneScan 500 LIZ size standard. Alleles were identified with GENEMAPPER v 4.0 (Applied Biosystems, USA).

| Population genetic diversity analysis
The number of alleles, observed heterozygosity (H O ), and polymorphism information content (PIC) were analyzed by the macros in Microsatellite Tools (Park, 2001). The null allele frequencies were estimated using the software FreeNA (Chapuis & Estoup, 2007).
Deviation from Hardy-Weinberg equilibrium (HWE) for each loci/ population combination, linkage disequilibrium (LD) among loci within each population, pairwise mean population differentiation (F ST ), and inbreeding coefficients (F IS ) were estimated in GENEPOP v4.0.11 (Rousset, 2008). We used the HP-RARE v1.1 (Kalinowski, 2005) to test allelic richness (A R ) and allelic richness of private alleles (P AR ) of each site. We used GENCLONE v2.0 (Sophie & Khalid, 2006) to estimate the total number of alleles (A T ) and the unbiased expected heterozygosity (H E ). We compared the number of alleles (A S ) and heterozygosity (H ES ) among samples with different sample sizes in GENCLONE.

| Population genetic structure analysis
First, population genetic structure was investigated using STRUCTURE v2.3.4 (Earl & vonHoldt, 2012). We used 30 replicates of each K value from 1 to 10, with 200,000 Markov chain Monte Carlo iterations and a burn-in of 100,000 iterations. The results were submitted to the online software Structure Harvester v0.6.94 (Earl & vonHoldt, 2012) to determine the optimal K value by a Delta K method. Membership coefficient matrices (Q-matrices) associated with the optimal K were processed using CLUMPP v1.12 (Jakobsson & Rosenberg, 2007), and then visualized using DISTRUCT v1.1 (Taubert et al., 2019).
Second, discriminant analysis of principal component (DAPC) was used to analyze population genetic structure under default settings (Jombart et al., 2008), complementing the STRUCTURE analysis. This analysis was run using an R package adegenet v1.4-2.

| Gene flow analysis
We used the BAYESASS v3.0.4 (Wilson & Rannala, 2003) to estimate the recent migration rates among populations of the two species.
First, we conducted preliminary runs (10,000,000 steps) to adjust mixing parameters for allele frequencies and inbreeding coefficients.
Then, we carried on ten longer runs of 100,000,000 steps using different start seeds with a sampling frequency of every 1000 steps.
The trace files of 10 longer runs were combined using Tracer v1.6 (Rambaut et al., 2018) to calculate mean migration with a burn-in of 50,000. We used the migration rates of the two weevils as measures of gene flow.

| Analysis of the influence of geographic and climatic factors on genetic variation
Geographic and climatic effects on population genetic variation were examined with two methods. First, the Mantel test was used to test the presence of isolation by distance (IBD) and isolation by

| Population genetic structure
Structure analysis showed that the optimal K was two for TTW and  (Figure 2b and c). Individuals from different geographic regions were genetically divided into clusters similar to those identified from the STRUCTURE analysis.
Overall, population genetic structure analysis showed that some populations from different geographic locations could be assigned to the same genetic cluster (e.g., BJHD, BJHR, NXLW, NXPL, NXZW, and SDTA in TTW, see Figure 1), while some populations collected from the same area could be assigned to different groups (e.g., BJHR and BJHD in TRW, see Figure 1). Additionally, individuals from the same population could be assigned to different clusters (e.g., TJTJ in TTW and NXZW in TRW).

| Dispersal of TTW and TRW in China
The distribution range of A. altissima matches the current distribution range of the two weevils, with all three taxa mainly distributed Ailanthus altissima may have originated from the Indian subcontinent, first spreading to Tibet with the collision of plates, and then spreading to other regions such as East Asia, Europe, and even North America via Tibet (Su et al., 2021). In addition, the humid and subhumid areas of the Qinghai-Tibet Plateau provide a suitable climate for A. altissima (Yan, 2019). Perhaps this climate region provides a source location for the two weevils, although samples were unavailable from this region.
In our study, TTW and TRW provide a pattern of cross-diffusion between the northeast and southwest regions. According to the 2014-2017 National Forestry Pest Survey Results (Song, 2020), TTW and TRW are distributed in northeast, northwest, and southeast China. However, we found that population diversity in some southeast distribution areas is lower than that of northwest and north China.
This might be caused by a lower population density in some southern populations which could account for sampling difficulties at some of our locations (e.g., we only collected five TTW individuals in HBJZ).

| The influence of climate and geographic factors on population genetic variation
The results of IBD and IBE suggest that the genetic distance of two species has no strong association with geographic distance or climate factors. Since the IBE analysis considered the correlation between two matrices using a Mantel test and reduced the dimensions of the original data, we conducted a multivariate RDA analysis. The results showed that the variation in genetic distance between samples may have some association to climate factors in both the TTW and TRW populations. Among these, temperature and precipitation may affect patterns of genetic variation in the two weevils. The importance of these variables has previously been noted for another pest, the melon thrips Thrips palmi (Cao et al., 2019).

| Implications for pest management and usage in biological control of invasive plant
The integration of DNA data and pest risk assessment offers important monitoring tools for evaluating the relative success of pest invasions (Stepien et al., 2005). Northwest China (including the Ningxia Hui Autonomous Region, Qinghai, Shanxi Province, etc.) has classified TTW and TRW as key prevention and quarantine targets in recent years. The host plant of TTW and TRW, the tree-of-heaven, is a plant that tolerates a wide range of climatic conditions and can be used when planted as seedlings to construct shelterbelts in northern and northwestern China (Hu, 1979). Our results point to high levels of differentiation among some populations but low genetic diversity, consistent with spread through human activities and a buildup of the pest from low initial numbers.
Due to the weak dispersal capacity of TTW and TRW, quarantine measures and disinfestation should be effective and used to limit transmission during seedling transfer between Chinese provinces and cities. Such measures should prevent the two species of weevil from spreading further.
In the United States, A. altissima has been listed as an invasive weed because it has destroyed the landscape of urban areas in recent years. The two weevils studied here are monophagous woodboring pests which only feed on A. altissima. The latitude and climatic conditions of the United States and China are similar, and there is a substantial area suitable for these two weevils in the United States (Ji et al., 2017). We propose that these two weevils could be spread through humans to act as biological control agents in controlling this tree.

| CON CLUS ION
In our study, we used microsatellite markers to explore the population genetic diversity and structure of two Eucryptorrhynchus weevils. We found that populations of both TTW and TRW have high levels of differentiation but low genetic diversity. The patterns of genetic structure and intraspecific admixture point to the human-mediated dispersal in these two weevils. We found that temperature and precipitation may indirectly associate with genetic differentiation in the two weevils. The results provide

CO N FLI C T O F I NTE R E S T
There are no conflicts of interest. F I G U R E 4 RDA analysis on genetic variance of Eucryptorrhynchus brandti (TTW) and E. scrobiculatus (TRW) explained by the environmental effects of climate and geography. A full RDA model was run by considering environmental and geographic effects simultaneously (a, TTW; c, TRW), and a partial RDA model was run by constraining geographic effects to analyze the correlation of environmental variables (b, TTW; d, TRW). Individuals from the same population are indicated by circles with the same color. PCNM1-2, geographic variables; the highest temperature of the warmest month (bio5), Min temperature of coldest month (bio6), temperature annual range (bio7), precipitation of wettest month (bio13), precipitation of driest month (bio14), precipitation seasonality (coefficient of variation) (bio15) and precipitation of wettest quarter (bio16), climate variables. Correlations of each variable are indicated by an arrow. Long arrows indicate a high correlation between the variable and genetic distance Writing -original draft (equal); Writing -review & editing (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
Microsatellite data generated in this study were deposited to the Dryad data repository (https://doi.org/10.5061/dryad.vx0k6 djtr).